better formulation of polyroot

05-open_channel.Rmd

@@ -392,10 +392,11 @@ For a 0.25 m rise, and using $q = Q/b = 2.2/0.5 = 4.4$, combining Equation \@ref
 $$E_2 = E_1 - 0.25 = 2 + \frac{4.4^2}{2(9.81)(2^2)} - 0.25 = 2.247 - 0.25 = 1.997 ~m$$
 From the specific energy diagram, for $E_2=1.997 ~ m$ a depth of about $y_2 \approx 1.6 ~ m$ would be expected, and the flow would continue in subcritical conditions. The value of $y_2$ can be calculated using Equation \@ref(eq:open-rect1):
 $$1.997 = y_2 + \frac{4.4^2}{2(9.81)(y_2^2)}$$ which can be rearranged to $$0.9967 - 1.997 y_2^2 + y_2^3= 0$$
-Solving a polynomial in R is straightforward using the `polyroot` function and using `Re` to extract the real portion of the solution.
+Solving a polynomial in R is straightforward using the `polyroot` function and using `Re` to extract the real portion of the solution (after filtering for non-imaginary solutions).
 
 ```{r openrect-poly, message=FALSE, warning=FALSE}
-Re(polyroot(c(0.9667, 0, -1.997, 1)))
+all_roots <- polyroot(c(0.9667, 0, -1.997, 1))
+Re(all_roots)[abs(Im(all_roots)) < 1e-6]
 ```
 
 The negative root is meaningless, the lower positive root is the supercritical depth for $E_2 = 1.997 ~ m$, and the larger positive root is the subcritical depth. Thus the correct solution is $y_2 = 1.64 ~ m$ when the channel bottom rises by 0.25 m.

docs/flow-in-open-channels.html

@@ -908,9 +908,10 @@ <h2><span class="header-section-number">5.6</span> Flow in Rectangular Channels<
 <span class="math display">\[E_2 = E_1 - 0.25 = 2 + \frac{4.4^2}{2(9.81)(2^2)} - 0.25 = 2.247 - 0.25 = 1.997 ~m\]</span>
 From the specific energy diagram, for <span class="math inline">\(E_2=1.997 ~ m\)</span> a depth of about <span class="math inline">\(y_2 \approx 1.6 ~ m\)</span> would be expected, and the flow would continue in subcritical conditions. The value of <span class="math inline">\(y_2\)</span> can be calculated using Equation <a href="flow-in-open-channels.html#eq:open-rect1">(5.20)</a>:
 <span class="math display">\[1.997 = y_2 + \frac{4.4^2}{2(9.81)(y_2^2)}\]</span> which can be rearranged to <span class="math display">\[0.9967 - 1.997 y_2^2 + y_2^3= 0\]</span>
-Solving a polynomial in R is straightforward using the <code>polyroot</code> function and using <code>Re</code> to extract the real portion of the solution.</p>
-<div class="sourceCode" id="cb43"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb43-1"><a href="flow-in-open-channels.html#cb43-1" tabindex="-1"></a><span class="fu">Re</span>(<span class="fu">polyroot</span>(<span class="fu">c</span>(<span class="fl">0.9667</span>, <span class="dv">0</span>, <span class="sc">-</span><span class="fl">1.997</span>, <span class="dv">1</span>)))</span>
-<span id="cb43-2"><a href="flow-in-open-channels.html#cb43-2" tabindex="-1"></a><span class="co">#&gt; [1]  0.9703764 -0.6090519  1.6356755</span></span></code></pre></div>
+Solving a polynomial in R is straightforward using the <code>polyroot</code> function and using <code>Re</code> to extract the real portion of the solution (after filtering for non-imaginary solutions).</p>
+<div class="sourceCode" id="cb43"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb43-1"><a href="flow-in-open-channels.html#cb43-1" tabindex="-1"></a>all_roots <span class="ot">&lt;-</span> <span class="fu">polyroot</span>(<span class="fu">c</span>(<span class="fl">0.9667</span>, <span class="dv">0</span>, <span class="sc">-</span><span class="fl">1.997</span>, <span class="dv">1</span>))</span>
+<span id="cb43-2"><a href="flow-in-open-channels.html#cb43-2" tabindex="-1"></a><span class="fu">Re</span>(all_roots)[<span class="fu">abs</span>(<span class="fu">Im</span>(all_roots)) <span class="sc">&lt;</span> <span class="fl">1e-6</span>]</span>
+<span id="cb43-3"><a href="flow-in-open-channels.html#cb43-3" tabindex="-1"></a><span class="co">#&gt; [1]  0.9703764 -0.6090519  1.6356755</span></span></code></pre></div>
 <p>The negative root is meaningless, the lower positive root is the supercritical depth for <span class="math inline">\(E_2 = 1.997 ~ m\)</span>, and the larger positive root is the subcritical depth. Thus the correct solution is <span class="math inline">\(y_2 = 1.64 ~ m\)</span> when the channel bottom rises by 0.25 m.</p>
 <p>A vertical line or other annotation can be added to the specific energy diagram to indicate <span class="math inline">\(E_2\)</span> using <em>ggplot2</em> with a command like <code>p1 + ggplot2::geom_vline(xintercept = 1.997, linetype=3)</code>. The <em>hydraulics</em> R package can also add lines to a specific energy diagram for up to two depths:</p>