I have the following ordinary differential equation for $y(x)\small x\in [0 \ h]$ where $h>0$
\begin{equation} Ty''-(P_{0}+\gamma x)[1+y'^{2}]^{\frac{3}{2}} = 0; \quad y(0) = 0;\quad \frac{1}{y'(0)}=0 \end{equation}
A related question can be found here, and subsequently this related question will called as the "original question" in this query. The solution to the ode is very similar to calculating the 2D shape of a liqid droplet resting on a solid surface due to surface tension. However, it comes from a geotechnical field problem & not fluid mechanics.
I have solved the above ode in MATLAB using the ode45 suite for given $T$, $\gamma$ and $P_{0}$. Now however, the objective has changed. One needs to find the the appropriate value of $T$, given a constraint that the length of the curve, $L$ plus the ordinate at the end of domain $[0 \ h]$ i.e. $L+y(h)$ is a given constant say $c>0$.
By the way, $h$ can be determined for real solutions by solving for the positive root of the quadratic equation $Ax^2+Bx-2=0$, where $A = \frac{\gamma}{2T}$ and $B=\frac{P_0}{T}$. Moreover, $y'(x) = 0$ at say $x=h_{0}$, provided a similar quadratic equation $Ax^2+Bx-1=0$ is solved for its positive roots. So you can see that given the the coefficients $T$, $\gamma$ and $P_{0}$, the solutions can be numerically evaluated easily.
However,now one needs to find the appropriate coefficients to satisfy a given constraint.
In other words: Given $P_{0}$ and $\gamma$, find appropriate $T$ such that $L+y(h)=c$ for the abovementioned IVP.
Before trying to find an optimal algorithm for calculating $T$, I want to calculate $T$ using a linear search. For this I seed $T$ with a guessed value. This guess is based on the problem's similarity with the case of a liquid droplet as mentioned earlier. For large values of $T$, I assume the shape will be circular. In that case $y(h)=0$ and $h=2r$, where $r$ would be the radius of the circular shape. Hence $h = \frac{2c}{\pi}$. Rearranging the first quadratic equation will solve for the guessed value of $T$, now that $h$ is known.
Now I start with this guessed value of $T$ and keep decreasing it at every iteration to end up with the appropriate value within a specified tolerance. However, there are issues
- Since ode has a singularity at $x=0$, I specify the IVP to start at $x=0.0001$ instead i.e. $y(0.0001)=0$. Does this detail effect the computation, I don't know. I doesn't seem to be in the original question
- Since rounding and precision are the evils of computation calculating & specifying exact $h$ is not possible. A bit more than the required $h$ , however, means complex solutions which is not what concerns me. So I specify while numerically integrating in ode45, the problem as the ode equation with domain as $[0.0001 \ h_{modified}]$ with $y(0.0001)=0$ as stated earlier and $h_{modified} = h-100\epsilon$. Note that, $\epsilon$ here is the machine precision. This procedure works for the original question. But for linear search, iterations invariably one end up in complex solutions or NaN.
- The iterations are supposed to be broken when an error term is less than a speciifed tolerance. This error is defined as $|c-(L+y(h_{modified}))|$ where $c$ is known while $L$ and $h_{modified}$ are computed quatities in each iteration.
Frankly, I want to know what am I doing wrong in the linear search that I end up getting complex solutions midway in the iterations. Is there a way to avoid or improve this approach/procedure.
Note 1: 0.0001 is an arbritraily chosen low value close to zero.
Note 2: $100\epsilon$ will guarantee that the rounding error is 100 units in the last place.
Note 3: Instead of choosing $100\epsilon$ one can choose some other value like $10000\epsilon$, but still it doesn't work in the linear search iterations.
Note 4: I computed $\epsilon$ from matlab function eps. Specifically $\epsilon = $eps(h) in matlab syntax.