# How to solve a nonlinear equation?

I'm going to ask my question, explain my problem, show how I thought to solve the problem and where I got stuck, and then re-state my question.

Question: how might I solve this equation:

$$\tag{1} Re_{p}(S)=\frac{S}{1+0.15Re_{p}^{0.687}}$$ for $Re_{p}(S)$?

Explanation of my problem: I am trying to follow a paper, it goes like this:

The settling velocity for a particle falling in a viscous fluid under gravity is $$\tag{2} v_g=\sqrt{\frac{4}{3} \frac{dg}{C_D} \left(\frac{\rho_p-\rho}{\rho}\right)}$$

where $C_D$ is the drag coefficient, a function of the particle Reynolds number

$$\tag{3} Re_p=\frac{v_g d}{\nu}$$

The correlation for $C_D$ is given as

$$\tag{4} C_D=\frac{24}{Re_p}\left(1+0.15Re_{p}^{0.687}\right)$$

The three equations (2)-(4) can be reduced to a single nonlinear equation with respect to the Reynolds number

$$\tag{5} Re_{p}(S)=\frac{S}{1+0.15Re_{p}^{0.687}}$$

where

$$\tag{6} S=\frac{d^3 g}{18 \nu^2}\left(\frac{\rho_p-\rho}{\rho}\right)$$

The authors of the paper then say the solution of Eqn (5) is shown in the Figure below (as the solid line), which allows one to find the settling velocity

$$\tag{7} v_g=\frac{\nu Re_p(S)}{d}$$

as a function of the settling number $S$.

How I thought to solve the problem and where I got stuck: I thought I might be able to solve for $Re_p(S)$ by rearranging Eqn (1) as so

$$\tag{8} S=Re_p(1+0.15Re_p^{0.687})$$

distributing

$$\tag{9}S=Re_p+0.15Re_p^{1.687}$$

and then procede by performing the "completing the square" method, e.g., dividing through by $0.15$:

$$\tag{10} Re_p^{1.687}+\frac{Re_p}{0.15}=\frac{S}{0.15}$$

and then completing the square by adding $\left(\frac{Re_p}{0.3}\right)^2$ to both sides:

$$\tag{11} Re_p^{1.687}+\frac{Re_p}{0.15}+\left(\frac{Re_p}{0.3}\right)^2=\frac{S}{0.15}+\left(\frac{Re_p}{0.3}\right)^2$$

However, since my first term on the LHS is not to the power of 2 I can't proceed to simplify the LHS.

Back to the question: How did the authors of the paper solve Eqn(1) to produce the solid line in the figure given above? (Solved for $Re_p(S)$)

(Also, please let me know if my question title is appropriate for the question I am asking)

• Are you sure they didn't solve it numerically? – stochastic Jan 5 '18 at 21:36
• There is no "completing the square" or related solution since the exponent is a numeric value. The solution will be done numerically. – Andreas Jan 5 '18 at 21:38
• @Andreas what sort of numerical method should I look up to begin to understand how I might solve Eqn (1)? – Armadillo Jan 5 '18 at 22:02

This is similar to solving equations like $\tan(x) = \log(x)$. There is no way to "isolate" for $x$. In order to solve for the Reynold's numer in your equation, an iterative (i.e. numerical) approach must be taken. I assume you are in engineering, so eventually you will take a Numerical Methods course where you will learn all about solving non-linear equations.

For now, you can take the following simple iterative approach. Recall your equation $$Re_p = \frac{S}{1 + 0.15Re_p^{0.687}}$$ Suppose you want to solve this for $S = 10^0 = 1$. Then, $$Re_p = \frac{1}{1 + 0.15Re_p^{0.687}}$$ Make an initial guess for $Re_p$, call it $Re_p^0 = 10$. Use the above equation to iteratively get better guesses for $Re_p$. That is, $$Re_p^{i + 1} = \frac{1}{1 + 0.15{Re_p^{i}}^{0.687}}$$ The following table summarizes the results for the first few iterations. $$\begin{array}{c | c} i & Re_p^i \\ \hline 0 & 10\\ 1 & 0.5772098\\ 2 & 0.9068355\\ 3 & 0.8770160\\ 4 & 0.8794782\\ 5 & 0.8792734\\ 6 & 0.8792904\\ 7 & 0.8792890\\ 8 & 0.8792891\\ 9 & 0.8792891\\ \end{array}$$ Of course, you'd want to automate this by writing a simple computer program.

As a final note, I suspect that the authors might have gotten values of $S$ by simply substituting values for $Re_p$ in the above equation, and plotting the result for a range of $Re_p$ that are of interest.

$\color{Red}{\textrm{Mathematica code:}}$

R[S_] := Rep /. First@NSolve[S == Rep + 15/100*Rep^(1687/1000), Rep, Reals];

r1[S_] := S; (*  a dummy function *)

LogLogPlot[{R[S], r1[S]}, {S, 10^-1, 10^3},
PlotRange -> {Automatic, {10^-1, 10^3/2}}, AxesLabel -> {"Rep", "S"},
PlotLabel -> "Settlement of spherical particle in viscous fluid",
PlotStyle -> {Blue, {Red, Dashed}}, PlotLegends -> {Automatic}]


$\color{Red}{\textrm{Maple code:}}$

restart:
R := proc (S) fsolve(S = Rep+0.15*Rep^1.687, Rep, maxsols = 1) end proc:
plots:-loglogplot(['R(S)', S], S = .1 .. 10^3, view = [0.1 .. 1000, 0.1 .. 500]);


$\color{Red}{\textrm{Maxima code:}}$

R(S):= find_root (S=Rep + 0.15*Rep^1.678, Rep, 0.10,999)$; plot2d ([R(S),S], [S, 0.11, 999],logx,logy)$


As already said, you need some numerical method to solve the nonlinear equation $$S=R+k\,R^{\,a}\qquad \text{with}\qquad k=0.15\qquad \text{and}\qquad a=1.687$$ and, as usual, if you want to save iterations, you need a "good" guess for starting.

Playing with inequalities, you know that the solution is somewhere between $$R_1=\frac S {1+k}\qquad \text{and}\qquad R_2=\left(\frac S {1+k}\right)^{1/a}$$

So, let us take as estimate the geometric mean of these bounds $$R_0=\sqrt{R_1 R_2}=\left(\frac S {1+k}\right)^{\frac{1+a}{2a}}$$

Now, consider that you look for the zero of the function $$f(R)=R+k\,R^{\,a}-S$$ $$f'(R)=1+a\, k\, R^{\,a-1}$$ Newton iterates will be given by $$R_{n+1}=R_n-\frac{R_n+k\, R_n^{\,a}-S}{1+a\, k\, R_n^{\,a-1}}$$

For illustration purposes, let us try with $S=10^3$; Newton method will provide the following iterates $$\left( \begin{array}{cc} n & x_n \\ 0 & 219.186 \\ 1 & 170.109 \\ 2 & 165.981 \\ 3 & 165.950 \end{array} \right)$$ which is the solution for six significant figures.

For sure, this was for computing a single value. If you want to build the curve, for a second point, use the previous result as estimate. Suppose that we want now $S=1100$; then, the iterates will be $$\left( \begin{array}{cc} n & x_n \\ 0 & 165.950 \\ 1 & 176.500 \\ 2 & 176.303 \end{array} \right)$$

For example, using $S=1$, Newton iterates would be $$\left( \begin{array}{cc} n & x_n \\ 0 & 0.894667 \\ 1 & 0.879285 \\ 2 & 0.879268 \end{array} \right)$$

Edit

You could do simpler (not to say much better) if you try to inverse the model. Generating $S$ from given $R$ and using least-square method to minimize the sum of squares of relative errors, with a very simplistic model, you could end with $$R\approx 1.35463\, S^{0.869369}-0.384677\, S$$ which would give pretty good guesses to be used for $R_0$.

For the worked examples, the values generated as estimates are $164.772$, $173.771$ and $0.969949$. Newton method will probably require a couple of iterations for convergence to high accuracy.

• This answer is incomplete as it lacks a proof of convergence of the Newton iteration, which can be easily remedied. – Hans Jan 7 '18 at 4:35
• @Hans. The first derivative is always positive; so only one root. – Claude Leibovici Jan 7 '18 at 4:46
• The uniqueness of the root of the original equation does not guarantee the convergence of the Newton iteration. – Hans Jan 7 '18 at 5:13